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Abstract 

Yukawa potentials are often used as effective potentials for systems as colloids, plasmas, etc. 
When the Debye screening length is large, the Yukawa potential tends to the non-screened 
Coulomb potential ; in this small screening limit, or Coulomb limit, the potential is long ranged. 
As it is well known in computer simulation, a simple truncation of the long ranged potential and 
the minimum image convention are insufficient to obtain accurate numerical data on systems. 
The Ewald method for bulk systems, i.e. with periodic boundary conditions in all three directions 
of the space, has already been derived for Yukawa potential [cf. Y., Rosenfeld, Mol. Phys., 88, 
1357, (1996) and C, Salin and J.-M., Caillol, J. Chem. Phys., 113, 10459, (2000)], but for 
systems with partial periodic boundary conditions, the Ewald sums have only recently been 
obtained [M., Mazars, J. Chem. Phys., 126, 056101 (2007)]. In this paper, we provide a closed 
derivation of the Ewald sums for Yukawa potentials in systems with periodic boundary conditions 
in only two directions and for any value of the Debye length. A special attention is paid to the 
Coulomb limit and its relation with the electroneutrality of systems. 
*Electronic mail: Martial.Mazars@th.u-psud.fr 
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1 Introduction 



Yukawa interaction potentials, named likewise after the achievement of the meson theory by 
Hideki Yukawa in 1935 [lj, were introduced by Debye and Hiickel in 1923 as a mean field 
approximation in the Poisson-Boltzmann equation for ionic solution [2j ; these potentials, in 
their static form, are given by 



where k is the inverse screening length, and Q an effective Yukawa charge. These quantities are 
related to some physical parameters of systems ; some examples are given in Table 1. 
The Yukawa potential is solution of the Helmholtz equation 



with the boundary condition V(oo) = 0. 

In theories of ionic systems as liquids [3j 2], colloids El [8] , plasmas [9j [10], etc., k^ 1 is 
known as the Debye length, this length gives an estimation of the radius of the screening sphere 
and Q is an effective charge that may be related to the surface potential ips of particles with 
hard-core interaction. More precisely, for hard spheres systems we have 



with a the radius of the hard spheres (cf. Table 1). 

In Yukawa's meson theory and according to the general principle of the quantum theory, k is 
related to the mass m n of the meson by m n = kU/c and Q is an effective coupling constant for 
nuclear interactions pQ. 

On general grounds, Yukawa interaction potentials can be considered as a reasonable approx- 
imation of effective interaction potentials between particles when some microscopic degrees of 
freedom may be approximated to a continuous background that screens the direct interaction 
between particles, while the spherical symmetry of the interaction is preserved [11] . Despite 
the apparent simplicity of the analytical form of the Yukawa potentials in Eq.(l), the physical 
mechanisms by which such effective potentials may be derived from the microscopic degree of 




(1) 




(2) 



Q = ips aexp(fca) 
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freedom can be quite complicated and these mechanisms depend strongly on the properties of 
studied system. For instance, one may compare the difference between the charging mechanisms 
of dust grains in dusty plasma systems [9l [TOT [T5| \W[ E] an d the equilibrium Poisson-Boltzmann 
approximation or the other approximations schemes in colloidal systems [UJ [TZj [TBI Q3] . 
In computer simulations, one uses periodic boundary conditions to avoid irrelevant surface bias. 
As the lattice sums of Yukawa potential are absolutely convergent for any value of the inverse 
Debye length k, one may safely use a truncation of the potential and minimum image convention, 
as long as k is large enough. But, in the Coulomb limit (k — ► 0), truncation of the potential 
introduces important bias and errors |18[ [T9l [20] ; therefore, there is a domain with fc/0 where 
one has to handle lattice sums of Yukawa potential with caution by using convenient algorithm 
as Ewald methods |21[ I22j. As outlined before by authors of refs.|21[ I22j. truncation of Yukawa 
potentials for bulk systems may safely be used as long as exp(— nL)/L is small. Typically, for 
systems with number density p* ~ 1 with N ~ 1000 Yukawa particles, truncation of the Yukawa 
potential can be used for k* = kg > 2 (cf. Table 1) |23 [ 1241 1251 26J. It is also worthwhile to note 
that a code with Ewald sums for Yukawa potential, for any value of n, can easily be obtained, 
with minimal modifications, from a code with Ewald sums for Coulomb potential |21[ I22j. 
Many interesting experimental systems in plasmas physics \27\ [28] and in colloids science are 
quasi- two dimensional systems |30[ [31] , or even quasi-one dimensional systems [32] . To study 
such systems with numerical simulations, one have to use partial periodic boundary conditions 
in one or two directions. The purpose of the present work is to provide a detailed derivation of 
the Ewald method for Yukawa potentials in quasi-two dimensional systems. It seems that Ewald 
method for Yukawa potentials in quasi two dimensional systems have already been used in some 
works of refs.|28|. as outlined in ref.[29] ; however, the analytical details of the Ewald sums for 
Yukawa potentials are not gven in these works. In a paper numbered II, the Lekner sums |33[I35] 
for Yukawa potentials will also be given and in the paper III, we give numerical results for a 
bilayer system of particles interacting via a Yukawa potential ; some preliminary results for this 
Yukawa bilayer system are given in the section 4 of the present paper. The analytical form of 
Ewald sums for Yukawa potential obtained in the present work shows that, as for bulk systems, 
a code can be easily written from one code for Coulomb potential with minimal modifications. 
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In this paper, a special attention is paid to the Coulomb limit and its relation with electroneu- 
trality. From this analysis we obtain the constants useful to compute the total energy. These 
constants depend on the geometry of systems and on the way the electroneutrality is achieved. 
The present work is organised as follow. In section 2, we derive the Ewald sums of Yukawa in- 
teraction potential for quasi-two dimensional systems by using a method similar to the one used 
by Grzybowski and co-workers for the Coulomb interaction potential [36]. At the end of section 

2, the particle-particle interaction energy for a system of N particles in a quasi-two dimensional 
system, this result does well agree with a recent result |37| obtained from both the Ewald sums 
for Yukawa potential in three dimensional systems [22] and the Parry method [38]. In section 

3, the small screening limit and its relations with electroneutrality of systems are outlined. In 
particular in this section, we show explicitly how singular contributions are cancelled by elec- 
troneutrality for three different systems ; these results allow to link the total energy of systems 
interacting with Yukawa potential to the energy of systems interacting with Coulomb potential. 
Section 4 is devoted to a discussion on the practical use of the Ewald method for Yukawa poten- 
tial and to numerical tests. For completeness two appendices have been added at the end of the 
paper. In appendix A, we provide a full derivation of the Ewald sums for Yukawa potentials in 
three dimensional systems with the method by Grzybowski and co-workers [36] ; the results of 
this appendix agree with the derivation done by Salin and Caillol who have computed the Ewald 
sums from the Green's function of the Helmholtz equation [22]. In appendix B, we compute the 
forces obtained with the Ewald method for molecular dynamics implementations. 

2 Ewald sums for Yukawa interaction in quasi-two 
dimensional systems. 

We consider a system made of N particles interacting via Yukawa potentials. The simulation box 
has partial periodic boundary conditions, with spatial periodicity L x and L y , applied respectively 
to directions x and y, whereas no periodic boundary conditions are taken in the third direction 
z, parallel to the unitary vector e z . In the simulation box, the position of the particle i is defined 



4 



by Ti and we set 

fi = 8{ -\- Zi e z (3) 
The particle-particle interaction energy is given by 

^ N N j N 

£ cc (Yukawa; re) = - E E QiQi^ij) + ~ E Q? $ o (4) 



1=1 r^j 1=1 



with $(r) defined by 



^ exp(-re | r + n |) ^ exp(-re | n |) 

n I + I n^o 1 n 1 

In previous equations, we have used the condensed notations 



| r + n |= y (x + n x L x ) 2 + (y + n y L y ) 2 + z 2 and | n |= \Jn 2 L 2 + 

with re x and n y integer numbers associated with periodic images of the box. In appendix A, one 

would find the computation for periodic images in all the three directions of the simulation box. 

The computation starts with the identity 

exp(-re|r + n|) 1 f°° dt re 2 2 

exp(- — - \ r + n\U) (6) 



| r + n | \/tt Jo \/i At 

that will allow to use later the Poisson-Jacobi identity. 
With Eq.(6), $(r) may be written as 

1 f°° dt , re 2 



*(r) =^=fJ=eM-^eM-\r + n^t) 
= ^?/J^ eXp( 4 )eXp( - |r + n|2t) 



(7) 



+ V^Jo ^M- Tt -zt)^eM-\s + n\ t) 



= Ii(r; a; re) + J 2 (r; a; re) 
where the integral has been split into a short ranged contribution I\{r; cc; re) and a long ranged 
contribution hir; a; re). In Eq.(7), a is the Ewald damping parameter and, for practical appli- 
cations, its value is chosen to balance efficiency and accuracy ; a comparison between the results 
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of appendix A and the results in ref.[22] shows that one may interpret the parameter a in Eq.(7) 
as the gaussian width of a charge distribution that screens the central charge in the short ranged 
contribution [Note : In ref.|22|. the parameter a is denoted by and At by a]. As remarkably 
outlined in ref.[22], for three dimensional systems, the optimum choice for a is independent of 
k and is the same as for Ewald sums for Coulomb interaction, i.e. as k — > 0. 
With the integral relation 

POC dt , k 2 n s r°° , , a 2 t 2 



^7T ex P(~77 -a 2 *) =1 dtexp(-^— - n 2 t 2 ) 
2 y/t 4t J2a 4 



(8) 



^— [ exp(Ka)erfc(aa H ) + exp(— Ka)erfc(aa )1 

2a 2a 2a 



the short ranged contribution to the potential is given by 



where 



1 ^ D(| r + n \;a; k) 



D(r; a; k) = exp(«;r)erfc(ar + — ) + exp(— «;r)erfc(ar — — ) (10) 

2a 2a 



On Figure 1, we show the function D(r; a; k) versus nr for several values of the ratio A = n/a ; 
on this figure, the function C(r; a; k), used to compute the short ranged contribution of forces 
and defined by Eq.(B.4) in appendix B, is also represented. 
The long ranged contribution is computed by using the Poisson-Jacobi identity 

£ exp(- | s + n | 2 t) = ±£ £ e lG - S «p(-^) (H) 
n G 

where A = L x L y is the surface of the periodic simulation cell ; thus l2(r;a;K) becomes 



r i ^ f a2 dx k 2 2 

Hr;a- K ) =—J o ^exp(---zx) 



> r »Gs a rf;c , G 2 + k 2 2 

+-V 2^ e / "T7^ ex P( -a z x ) 

A Zr 1 Jo x 3 / 2 4x ' 

G^o 

With the help of Eq.(8) and after some algebra, we obtain 



(12) 



/ -^rexp(— z 2 x) = [exp(/cz)erfc(— + az) + exp(— &z)erfc( az 

Jo x ' 4x k 2a 2a 

/7F F(k; z; a) 
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(13) 



and thus 

J 2 (r; a; k) = - [f(k; z; a) + £ e iG -*F(VG 2 + k 2 ; z; a)} (14) 
A G^o 

Therefore, $(r) is given by 

= 1 2 D(I ; + ^ I; " ;K) + ~ A [F(«; z; a) + £ F( VG 2 ^; *; a)] (15) 

2^ |r + n| 7l L ^ 

with D(| r + n |; a; k) and F(k; z; a) defined respectively by Eqs.(lO) and (13). 
The lattice sums in Eq.(5) are absolutely convergent, as a consequence, as long as k / 0, 
there is no diverging contribution in the analytical decomposition of <E>(r) into short and long 
ranged contributions. As k — > 0, a 1/K-asymptotic behaviour is found ; this stems from the 
conditionally convergent property of the Coulomb lattice sums. This diverging behaviour of the 
Coulomb lattice sums is handled by taking into account the boundary conditions as | r \— > oo 
and the electroneutrality of the simulation cell [36J . The analysis of the asymptotic behaviour of 
Yukawa-Ewald sums as n — > and their links with the Coulomb-Ewald sums will be discussed 
extensively in the next section. 

To compute the energy E cc of the interaction between the particles, the self contribution $o is 

needed. This contribution is computed as for <3?(r). We have 

exp(-K | n |) 1 ^ f°° dt k 2 2 

*o = E i n = 7^ E / 2 7 =exp(--)exp(- \ n\ 2 t) 

(16) 

H — = / —p exp( — z t) > exp(— n i) = / — = exp( — z t) 

where the last integral has been added and substracted to the long ranged contribution to permit 
the use of the Poisson-Jacobi identity. Then 

*o = I? ] +4 0) (17) 

and after simple algebra we find 

(0) 1 ^ D(| n \;a;K) 



2 r£o l n 



(0) _ 2tt ^ erfc(/G 2 + k 2 /2a) r a k 2 » 

Gr 



(18) 



With Eqs.(4), (15) and (18), the particle-particle interaction energy is given by 
£ cc (Yukawa; = ~ E QiQj E ~ ^ + " |; ° ;i6j + £ Q iQj £ F(v /^T^; 2 

iV AT 2 

+^4 E QiQjF(n; zy; a) - ( E 0?) (7= ex P("^2 ) " ^erfc( K /2a)) 

i,j i=l V 

(19) 

where the prime in the summation over the image indicates that the term i = j has to be omitted 
for n = 0. 

The total particle-particle interaction energy is given by E cc (Yukawa; k) and may be used in 
Monte-Carlo simulations ; Eq.(19) agrees with the results obtained by using another method in 
ref.[37]. For Molecular Dynamics simulations, the computation of forces is needed, this deriva- 
tion is done in appendix B from the results of section 2 and 3. 

For repulsive Yukawa interaction between particles, one needs to compute the particle-background 
and the background-background energies to obtain the total energy of the system and thus to 
be able to implement correctly a Monte-Carlo sampling of the phase space. These computations 
are done in the next section for three different neutralizing backgrounds, and their limits as 
k — ► to a non-screened Coulomb system are outlined. 



3 The Coulomb limit and electroneutrality. 

For Coulomb interactions in systems with periodic boundary conditions, the electroneutrality 
property of the system is important to allow to define the electrostatic potential. Electrostatic 
Yukawa interactions are equivalent to non-screened Coulomb interactions as the screening length 
tends to infinity (or the inverse screening length k — > 0) ; therefore, in the Coulomb limit, 
the Yukawa interaction energy of the system, computed with Eqs.(4-5), has to be equal to 
the Coulomb interaction energy. The other limit, k — > 00, is either a hard sphere system 
(a, a, Oij 7^ 0) or an ideal gas limit (a, a, Oij = 0). 

The series in Eq.(5) are absolutely convergent when k > (Yukawa) and conditionally convergent 
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when k = (Coulomb). For Coulomb interactions, the conditional convergence property of the 
series leads to a diverging behaviour and this diverging behaviour is cancelled by using the 
electroneutrality of the system ; thus some constants, that depend on the geometry of the 
simulation box and on the manner in which electroneutrality is achieved, arise. As long as 
k ^ 0, the series are absolutely convergent and no diverging behaviour is found in the total 
interaction energy. To achieve consistency between Yukawa and Coulomb potentials in the 
Coulomb limit, one has to recover from Eqs.(19) and (5) a diverging behaviour corresponding 
to Coulomb potential and one has to be able to cancel it by using the electroneutrality of the 
system. This section is devoted to the examination of the Ewald sums consistency between 
Yukawa and Coulomb potentials for quasi- two dimensional systems. 

In subsection 3.1, the diverging terms appearing in the particle-particle interactions (Eq.19) as 
k — ► are displayed and in subsections 3.2 to 3.4, we provide the results obtained for several 
backgrounds necessary to fulfill electroneutrality : a monolayer in subsection 3.2 (Fig 2. a), a 
bilayer (subsection 3.3, Fig 2.b) and a slab (subsection 3.4, Fig 2.c). For any other neutralizing 
backgrounds, one may also compute the particle-background and the background-background 
energies according to similar methods. 



3.1 The Particle-Particle interaction energy in the Coulomb 



According to Table 1, if temperatures are high and/or if counterions are low, the inverse Debye 
length may become quite small. In the special case k = 0, one should recover the properties 
of One Component Plasma systems, thus, as k — > 0, Coulomb interaction energy has to be 
recovered from Eq.(19). 
As k — > 0, we have 



limit. 



D(| r \;a;n 



) =2 erfc(a | »* |) + |erfc(a | r \) 



exp(— a 2 | r | 2 ) 



r | 2 k 2 + o(k 4 ) 




< 



(20) 



F(k; z-a) =2-2(2 erf(az) + 



1 



K + K 2 + o(K 3 ) 




e 
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we have also 

G G 

G F(V G 2 + k 2 ; z;a) = \ exp(Gz)erfc( h az) + exp(— Gz)erfc( az) 

2a 2a 

k 2 G G 2G G 2 

+^7 [(1 " Gz)e Gz erfc(£- + az) - (1 + Gz)e~ Gz erfc(^ - az) - — exp(--^ - a 2 z 2 )} 
2G 2a 2a a^/n 4a z 



and 

a k 2 , n „ . k , a k 3k 2 



(21) 



exp(--^r) - - erfc( — ) = - - + V + o(k 3 ) (22) 



Therefore, we have 



£ cc (Yukawa; k -»■ 0) = ^(Coulomb) + ^( £ Q^) 2 - + ( ^ Q 2 )^ + ( K 2 ) (23) 



with 



s cc (couiomb) = if; q % q 3 y: / erfc { a r |r ;^ n|) + £t £ g,g, £ ejG?Sl3 F ( G ; «) 

N 1 2 2 1 

~\ £ Qi^i [*j erf(a^) + — =e~ Q2 4J _ " ( £ Q 2 ) 

(24) 

the Coulomb particle-particle interaction energy for a system with quasi-two dimensional geom- 
etry [36]. In Eq.(23), the singular contribution as 1/k appears explicitly. This singular term is 
cancelled if the system of particles has the electroneutrality property (i.e. 2~2i Qi = 0); then the 
small screening limit (k — ► 0) of the energy of the system, i? cc (Yukawa; k) in Eq.(19), is exactly 
E cc (Coulomb) in Eq.(24), the energy of the system of particles in interaction via a Coulomb 
potential. The electroneutrality implies that some pairs of particles in the system have attrac- 
tive Yukawa interaction as, for instance, in the restricted primitive model of electrolytes (RPM) 
with Coulomb poyentials. 

However, for applications to plasmas or colloids, one rather uses the Yukawa One Component 
Plasma model (YOCP) in which all particles have effective charge of the same sign and interac- 
tions between particles are purely repulsive ; therefore the electroneutrality can not be fulfilled 
by summation over the particles. For YOCP with periodic boundary conditions in the three 
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directions of the space, the electroneutrality is obtained by embedding the particles in a uni- 
form neutralizing background [22] (cf. Appendix A). For quasi-two dimensional systems, such 
uniform volume neutralizing background can not be used since boundary conditions in the z- 
direction depend on the model studied. In the next three subsections, the small screening limits 
of three different neutralizing backgrounds, adapted to the quasi-two dimensional geometry, are 
examined. These backgrounds are monolayer, bilayer and slab neutralizing backgrounds, they 
are schematically represented on Figure 2. 

3.2 The monolayer neutralizing background 

The monolayer neutralizing background represented on Figure 2 (a) is a plan with an uniform 
surface charge density localised at z = 0, particles may be localised on both sides of the plan or 
their location may be restricted also to only one half-space. Hereafter, system (a) designs the 
system made by the monolayer and particles ; for this system, the charge density in the right 
hand side of the Helmholtz equation, Eq.(2), is given by 




(25) 



where 8{z) is the Dirac distribution. Then, the electroneutrality reads as 



NQ + Aa = 



(26) 



The energy of the system (a) is given by 



E^ a > (Yukawa; k) = i? cc (Yukawa; «) + E cB {n) + E B b(k) 



(27) 



with 




NQ 



A 




exp(— k | T*i — r + n |) 



r i — r + n 



< 



(28) 




N 2 Q 2 r 
2A 2 Js 




exp(— k | r' — r + n |) 



r' — r + n 
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E c b(k) is the interaction energy of the particles with the background, while Ebb(k) is the self 

energy of the background. 

After elementary algebra, we find 

(29) 

N 2 Q 2 r- f°° dt 



N*Q< r f°° dt k\ 



and thus we have 

NO 2 1 N N 2 2 1 

£ (a) (Yukawa; k) = £ cc (Yukawa; k) - 2tt— f — ^ exp(-K | Zi \) + vr — j- — (30) 

J\ Hi . ^ J\ Hi 

1=1 

In the Coulomb limit, all singular contributions as 1/n are cancelled by taking into account the 
electroneutrality (cf.Eq.(26)) ; in this limit we have 

NO 2 N 

£ (a) (Yukawa; k -> 0) = ^(Coulomb) + 2tt—j- J2 I *i I (31) 

A i=i 

The second contribution in Eq.(31) stems from E c b(k — > 0). The energy of the Coulomb system 
is obtained by taking k = 0. The energy i?( a ) (Yukawa; k) computed with the Ewald sums is 
therefore defined for all value of k and the limit as k — > is well defined : it corresponds to a 
quasi-two dimensional coulombic system made of N particles carrying charge Q (One Component 
Plasma), and a monolayer with a surface charge density a to achieve electroneutrality, Eqs.(25) 
and (26). 



3.3 The bilayer neutralizing background 

The bilayer neutralizing background is represented on Figure 2 (b), this background is made of 
two plans, h is the separation between both plans and a the surface charge density in each plan. 
Particles may be localised between the plans or outside the slab. We design by system (b) the 
system made by the bilayer and the particles. For this system, the charge density appearing in 
the Helmholtz equation is 

Pb(r) = J2 QA r ~ r i) + ~ \) + *8(z + \) (32) 

i 
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and the electroneutrality reads as 

NQ + 2Acr = 

For system (b), the energy is given by 

E^ (Yukawa; re) = ^(Yukawa; re) + E c b{k) + Ebb(k) 

with the same definitions as for system (a). With simple computations, one finds 

tt , ^ NQ*("f°°dt h . 2 re 2 . 

Ecb(k) = — y o ^372 exp(-(^ - - ) t - — ) 



+ 



=1 io ^ eXp( - (Zl+ 2 )t -^ 



^(«) ^{2 1 ^exp(- ¥ ) + y o _exp(-^-- 



4t 



then 



A^Q 2 1 N 11 
E^ (Yukawa; re) = S cc (Yukawa; re) — ir — ^ ( exp(— re | Zi — — |) + exp(— re | Zi + — |)) 

2 = 1 



(33) 



(34) 



(35) 



h 



N 2 Q 2 l f , ,,x 

+vr^--(l + exp(-re/i)J 



(36) 



Similarly to system (a), the singular contributions in system (b) are cancelled by the electroneu 
trality. As h — > 0, we recover 

E {b) (Yukawa; re) -► £ (a) (Yukawa; re) 
and as re — > and h 0, one has 



(37) 



(Yukawa; re ->■ 0) = ^(Coulomb) +tt^- ^(| + - | + | |)-7r^^-/i+o(re) (38) 



2 N 



h 



2/n2 



i=i 



i4 



The same analysis may be done easily for multilayer neutralizing backgrounds. 

If Nq = N/2 particles are confined in each plan, then we have Z{ = —h/2 or z,- L = h/2 and the 

energy of the system is 



(Yukawa; re -► 0) = £ cc (Coulomb) + 2ttLV/j + o(re) 



(39) 
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then E^ b > (Yukawa; k —* 0) is exactly the energy of the bilayer Wigner system studied in refs. 
I39j . This will provide an useful reference point. 

3.4 The slab neutralizing background 

The slab neutralizing background is represented on Figure 2 (c), the slab between the plans 
located at z = —h/2 and z = h/2 is filled with an uniform volume charge density p$. Outside 
the slab, the volume charge density is null. Particles may be localised inside or outside the slab. 
We design by system (c) the system made by the slab and the particles. For this system, the 
charge density appearing in the Helmholtz equation is 

Pc(r) = £ QiS(r - n) + P0 (8(z + |) - 9(z - £)) (40) 

i 

with 6{z) the Heaviside distribution. For systems (c), electroneutrality reads as 

NQ + Ahp = (41) 
Similarly to systems (a) and (b), the total energy of system (c) is given by 

£ (c) (Yukawa; k) = £ cc (Yukawa; k) + E cB (k) + E bb (k) (42) 



with 



ei / \ NQ 2 /•! r°° dt . , .o k\ 

E cB (.) =-—yftY,]_^dz] Q wmp (-(zt-z)t--) 

(43) 

E BB ( K ) = J g d8' J\ dz'(l - e-^cosh^')) 

To compute E cB (k), one has to consider the set of particles that are confined in the slab and 
the set of particles being outside the slab, therefore we define 

I = { particles i :| z% |> |} and J = { particles j :| zj |< |} 

and we note 

Card/ = iVj and CardJ = Nj 
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obviously we have Nj + Nj = N. With these definitions of sets I and J, we find 



E c b(k) 



NQ 2 v 4tt 
Ah K 2 



sinh(^) e~ K[z ' 1 +Nj- e~ Kh/2 £ cosh^) 



(44) 



_ . , N 2 Q 2 2vrr 2 _ Kh/2 . nh. 

Again, as n — > 0, the singular contribution is cancelled by the electroneutrality and we have 
(Yukawa; k > 0) = ^(Coulomb) + vr^^ _ I + J_ ^ | ^ | ^ ^] + 



2L 2 L iV 3 Nh 

As /i — > and with the constraint poh = a, we find easily 

£( c ) (Yukawa; «) -> (Yukawa; as) 



(45) 



(46) 



4 Discussion 

A comparison between Eq.(19), that gives the total particle-particle interaction energy for 
Yukawa potential, and Eq.(24) giving the energy for Coulomb potential, shows that a code 
for Yukawa potential in quasi-two dimensional systems may easily be built from existing com- 
puter codes for Coulomb potential. Apart from constants, the main modifications to perform 
are : to change the function erfc(a | r^- + n |) in the short range contribution to the Coulomb 
energy by the function D(| + n \;a;K)/2 ; in the long ranged contribution, for G ^ 0, to 
change G in the function F(G; z^; a) by VG 2 + k 2 and to modify the contribution for G = 0. 
When all constants arising from the backgrounds contribution and self energies are included, 
then we may obtain a code for Yukawa potential for any value of the inverse screening length 
K, and this code will be consistent with the Coulomb limit. One should note also that the 
complexity of the long ranged contribution is the same as for Coulomb potential : the sum over 
the pairs and reciprocal lattice vectors G may not be simplified into summations over particles 
(with exceptions for some particular systems as the bilayers studied in refs.[39l [35]). 
As already noted by authors of refs. [21l 122} I24j. the use of Ewald sums for Yukawa potentials is 
only necessary when k is small, otherwise a brute truncation of the potential and the minimum 
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image convention are sufficient to perform accurate computations. To be more specific, we may 
estimate a lower bound n cu t for k that allows to obtain a relative accuracy of the order 5 in the 
computation of the energy. The criterion we use to estimate K cu t is as follow. 
Consider a system of Yukawa particles with hard sphere diameter a, then the interaction energy 
at contact is E c = Q 2 exp(— na)/a, while, if the separation between particles is L/2, with L 
the smallest spatial periodicity of the simulation box, the energy is El = 2Q 2 exp(— kL/2)/L. 
Then, a direct truncation of the potential will allow to obtain a relative accuracy of the order 5 if 
El/E c < 5. Then, we obtain an estimation of the lower bound n cu t that allows safe truncations 
of Yukawa potentials for distance greater than L/2, as 

K * = Ka > K* cut = -1— ln(-|-) = ln(^|) (47) 

where we use the reduced units L* = L/a, k* = Ka, a* N = a 2 <JN and <jn is the surface coverage 
of the simulation box defined by <tat = N/L 2 , N being the number of Yukawa particles. 
On table 2, we give some value of K* cut for different values of the surface coverage and the number 
of particles, for relative accuracy 5 = 10 -4 and 10~ 6 . 

In Table 3, we report some numerical tests of the Ewald sums for Yukawa potentials. For these 
tests, we consider a pair of Yukawa particles with effective charge Q\ = +1 and Q% = —1 in 
a box with L x = L y = 1 (in these reduced units, we have kL = k = k*L*) ; for this system, 
one has Q\ + Qi = 0, thus we can use straightforwardly Eq.(19) to compute the energy of the 
system and as k — > 0, Eq.(23) is well behaved. The five configurations reported in Table 3 have 
already been considered before to test Lekner sums for Coulomb potentials [341 [33"1 [4^1 14T] ; the 
values reported in the Column 'Coulomb Limit 7 have been computed with analytical, Lekner [33] 
and/or the Hautman-Klein methods \4:2\ 140] . For k = 10 -6 , the results obtained with Eq.(19) 
perfectly agree with the Coulomb limit for Ewald damping parameters a = 6 and a = 12 and 
with k x k = 16 x 16 for the number of vectors G in the reciprocal space. The results obtained 
with a direct truncation of the Yukawa potential do not allow to reproduce the correct Coulomb 
limit and energies computed with a direct truncation for n = 1 are very different than the value 
obtained with Ewald sums. On Figure 3, we represent, for three configurations of the pair, the 
particle-particle energy E^Yukawa; k) as a function of k (symbols) and the energy computed 
with a direct truncation of the Yukawa potential (lines) ; the Coulomb limit for each configu- 
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ration is represented by a thick red line. It clearly appears on figure 3, that if k < 10 the use 
of Ewald sums is necessary ; this agrees with the results of Table 2 (k = k*L*). One may also 
note on Fig. 3 that, for the configuration (0.1,0.1,0.1) of the pair, by chance the energy computed 
with Ewald sums agrees quite well with the one computation using a direct truncation of the 
potential ; however an inspection of the numerical values shows that the difference is greater 
than 10~ 4 . 

On Figure 4, we represent AE cc (k, a) = E cc (K,a) — E cc (n,a = 6.) as a function of the Ewald 
damping parameter a for the three configurations of Figure 3. The contributions of the recipro- 
cal space are computed with k x k = 16 x 16 for the number of vectors G. These data show the 
stability of the Ewald sums for Yukawa potentials versus the arbitrary choice of the damping 
parameter a ; more precisely, an inspection of the numerical data shows that, if 4. < a < 15., 
then | AE cc (k, a) |< 10~ 6 . Therefore, a choice aL ~ 5. — 6. is convenient to use Ewald sums 
for Yukawa potentials, as already noted for Coulomb potentials |22j . On the one hand, if a is 
too small, then the minimum image convention is insufficient and additional images have to be 
included in the real space contribution ; on the other hand, if a is too large additional G-vectors 
have to be included in the reciprocal space contributions (cf. Fig. 4). 

On Table 4 and Figure 5, we report some preliminary results from Monte Carlo computations 
on Yukawa Bilayer systems. This system is the same as the one studied in refs.[35j [39] but 
the particles interact with a repulsive Yukawa potential and the electroneutrality is achieved 
with the background of subsection 3.3. The geometrical and physical parameters of the bilayer 
correspond to run a of ref.[35] ; more precisely, we have h = 1, N = 512, L = 28.36 and 
the coupling constant T = Q 2 jkTa = 196. Intralayer and interlayer energies and correlation 
functions are defined as in ref.[35]. On Table 4, data for Coulomb interaction are extracted 
from ref.[35], the results for Yukawa interaction have been obtained with the Ewald sums given 
by Eqs.(19) and (26). For k = 10 -4 , energies and correlation functions, on Fig. 5 (a,b), agree 
perfectly ; this shows that the Coulomb limit is correctly obtained with the results of section 2 
and 3. For larger values of k, the results obtained with the Ewald sums are compared to results 
obtained with a direct truncation of the Yukawa potential. For k = 0.1, energies and correlation 
functions obtained with a direct truncation differ significantly from the results obtained with 
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Ewald sums. For k = 0.25, energies obtained with a direct truncation of the Yukawa potential 
agree with energies computed with Ewald sums, but, as it appears on Fig. 5 (c,d), correlation 
functions disagree. For larger values of k, both energies and correlation functions computed 
with Ewald sums or with a direct truncation of the interaction potential are in agreement with 
a good accuracy. More details and properties on Yukawa bilayer systems will be given in the 
forthcomming paper numbered III. 

There are some deviations from Yukawa potentials [3j El H3] ; this is particularly true for quasi- 
two dimensional systems since counterions and coions density profiles in the z-directions are 
non-uniform. However, if, as in the theoretical analysis done in ref.[T4], one defines an effec- 
tive screening length that depends only on the z-coordinates of the pair of particles, then one 
may use the Yukawa- Ewald sums, Eq.(19), in replacing the uniform screening parameter k by 
a non-uniform k(z,z') (cf. Eqs. (25,26) of ref.[14j). Otherwise, if the screening of the Coulomb 
interaction is more complicated, one should come back to the reliability of the Yukawa interac- 
tion potential for the system considered. 

As a final comment to this paper, we would like to outline the fact that Coulomb limits may be 
obtained accurately for Yukawa potentials with the use of Ewald sums ; but, as it is shown in Sec- 
tion 3, the effective Coulomb limit of a system is depending on the manner the electroneutrality 
is restored from the background. 
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Appendix A: Ewald sums for Yukawa potential in three 

dimensional systems. 



To compute the particle-particle interaction energy of the system with periodic bound- 
ary conditions in three dimensions, we follow the same method as in section 2 for quasi-two- 
dimensional systems. Potentials are still defined by 

_ , . exp(-K I r + n I) . ^ ^ exp(-K | n |) 

*( p ) = E i r l n i and *° = E \ n ( ) 

but summations over the periodic images are taken along the three dimensions of the space. With 
periodic boundary conditions in all three directions of the space, the Poisson-Jacobi identity 
reads as 

E ^ V+n?t = Ujf /2 E ^ -P (-~) (A.2) 
n G 

From Eqs.(A.l) and (6), we have 

*W =^E/J^-p(-^)exp(-|r + n|^) 

^-P(-^-A)E-P(-|r + «| 2 t ) (A.3) 

= /{ 3D) (r; q; k) + 4 3D) (r ; q; k) 
and with Eq.(8), we find easily 

r(r ! °i.)^£ D(l ; r + ;f° i ' ) 

with D(| r + n |;a; k) given by Eq.(lO). I^^if] ct; k) is computed by using the Poisson-Jacobi 
Identity (A.2), we have then 

*«w - ^ ♦ 1 e - ( ^r 2) (-> 
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Without any difficulties, $o is also computed as follows 

1 ^ D(| n |;g;/c) 4vr ^ exp ( - (G 2 + K 2 )/4a 2 ) 4vr e"^/ 40 



(A.6) 



-2[4=exp(— ^)--erfc(— )1 
Therefore the particle-particle interaction energy is given by 
^-'(Yukawa; «) = \ £ Q, Qj VJ - gll^L^ 



(G 2 +K 2 )/4a 2 JV 2p -« 2 /4a 2 

*J G^o * =1 



AT 2 

EQ 2 )[^^p(-^)-i^c(^ 

„2 



The singular term as k — > in Eq.(A.7) is as l//e , in the Coulomb limit we have 

f(X» 2 ;HE<5?)f 



£<f '(Yukawa; K - 0) = ^'(Coulomb) + ^( £ + ( £ (J?) | + „(k 2 ) (A.S 



with -Ecc D ^ (Coulomb) given by 

^•(coulomb) ^^ms ^i';. 1 !-^"' ^!:^ £ e iGr « ^ 



£a 2 ) 



-7=' 

(A.9) 

The singular term in Eq.(A.8) has to be cancelled by the electroneutrality of the system. As- 
suming that for all particles we have Qi = Q > and that the electroneutrality is achieved by 
an uniform background po, then the charge density in the simulation box is 

P(3D)(r) = £Qi<S(r-r i ) + / o ( A - 10 ) 

i 

and the electroneutrality reads as 

NQ + Vp o = (A.ll) 
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Thus, for the three dimensional system, the total energy is 

(Yukawa; k) = E^P^ (Yukawa; k) + E% D) (k) + E$jj\n) 



with 



Ecb(k) = ~ 



NQ 
V 



2 N 

i=i JV n 



exp(— k | Ti — r + n 



ri — r + n 



exp(— k | r' — r + n 



2V Jv Jv ^ | r' — r + n 
After a simple computation of gaussian integrals, we obtain 

NQ 2 1 



Ecb(k) 



-Ait 



E BB U) =2tt 



V K 2 

NQ 2 1 



V K 2 

and therefore, for a system with an uniform neutralizing background, we have 



£ (3D) (Yukawa; k -> 0) = £ cc (Coulomb) + o(k) 
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Appendix B: Computation of forces. 



For molecular dynamics implementations, the computation of forces is needed ; in the following 
we give the forces obtained for Yukawa potentials in quasi-two dimensional systems. 
The force acting on particle i is given by 

Ft = if c) + Ff ] = -VrM CC) ~ VrA CS) 

(B.l) 



ET7i(cC; short) F (CC; long) jp(cS) 



where and E^ cS ^ are respectively interaction energies of particle i with other particles (cf. 
Eq.(19)) and with the neutralizing background (cf. Eqs. (30,36) or (43)). In Eq.(B.l), the force 
acting on particle i due to other particles is split into a short and a long ranged contribution. 
With the chain rule 



df(\ rjj + n |) _ d\r ij + n\ df{r) 

X 



(B.2) 

r=\r ij +n\ 



dxi dxi dr 

and Eq.(19), the short range contribution of the force acting on particle i due to particle j is 
given by 

F(cC; short) r-r t-i(cC; short) 

. /• — — V f ill ■ j- 

j/l l 3 /l 

4 n \ r ij + n \ 

(B.3) 

with Xij = Xi — Xj and F(r; a; k) defined by 

C(r; a; k) = (1 — nr) exp(rtr)erfc(ar + — ) + (1 + kv) exp(— Kr)erfc(ar — — ) (B.4) 

2a 2a 

Similarly, the long range contribution to particle-particle forces is given by 



jp(cc; long) _ 77i(cC; long) _ Z?( CC; lon s) _ ^ Tp( cc '< lon g) p, 

j/i ~ T i j/i ~ S i j/i Q z . j/i 



(B.5) 



From Eq.(19), we find 



- v Sl i-r: 1 ; lons) = -^jQiQj E MG-sa) f(V&T^-, Zij - a) g (b.6) 

G^o 
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and 



d_ E ^ long ) = _1_q.q. £ e iG. Sij B(v ^T^; Zij ; a) (B.7) 

with 



A; fc 

B(k; z; a) = exp(/cz)erfc(- h az) — exp(— fcz)erfc(- az) (B.8) 

2a 2a 

The force induced by the background on particle i depends on the geometrical parameter of 
the background. For systems (a), (b) and (c) described in section 3, there are some discontinu- 
ities in spatial distributions of the neutralizing background, hence forces are also discontinuous 
at those locations. Such behaviours may induce some complicated bias in molecular dynamics 
computations, especially when particles cross the surfaces of discontinuity. To implement molec- 
ular dynamics codes for systems (a), (b) and (c), one would have to estimate numerically the 
influence of these discontinuities on the trajectories of particles. 
Forces induced by backgrounds on particle i may be easily computed with 

F f) = -Lrt<*> ... e z (B.9) 

(a,b,c)/i Q z . (a,b,c)/i 2 v > 

For system (a), in the half-space Zj > 0, we have 

F\f /t = - 2 -^eM-^i)e z (B.10) 
for system (b), in the slab —h/2 < z; t < h/2, 

*S = sinh{KZi) e z (B.ll) 



(b)/i ~ A 

and for system (c), in the slab —h/2 < Zi < h/2 



I ON NO 2 p~ Kh l 2 

F t)l = ^^r^h- siQh ^ (B.12) 
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List of Tables 



Table 1 : Definitions of the Yukawa parameters Q and k in Eqs.(l,2) as functions 
of physical parameters and thermodynamical states of systems with screened- Coulomb 
interactions. For all systems, ks is the Boltzmann constant, T the temperature, e the 
elementary charge, eo the vacuum dielectric constant and e the dielectric permittivity of 
the medium. In the Debye-Hiickel approximation of electrolytes, qi and pi are respectively 
the charge and the number density of ions of species i, and for the Primitive Model of 
electrolytes, Oij = [pi + <Tj)/2, is the mean ionic diameter of ions i and j [H 0]. For 
colloidal systems, Z is the charge number and a the hard-core diameter of macroions, p 
is the number density of the monovalent counterions and A# = ^jek^T is the Bjerrum 
length [H [6], [3, E] ■ In Plasmas and Dusty Plasmas, the parameters for screening-Coulomb 
potentials are Z and a the charge number and the radius of the impurity (dust particle), 
n Q is the plasma density, Tj and T e are the ion and electron temperatures P [TO] • 

Table 2 : Values of n* cut for several values of the surface coverage and number of 
particles in the simulation box. If k* — Ka > K* ut , then one may use a direct truncation 
of Yukawa potential for distances greater than L/2 ; otherwise the long range of the po- 
tential has to be handled with a convenient algorithm. 

Table 3 : Numerical tests of the Coulomb limit of the Ewald sums for the Yukawa 
potentials. The configuration of the pair of particles is defined by (x\2, yi2, £12) particle 1 
with a charge Qi = +1 located at (0,0,0) and particle 2 with a charge Q2 = —1 located at 
(^12) 2/12? -212) in a box with L x = L y = 1. The values of particle-particle Coulomb energies 
are extracted from original works by others, in the original works the computations of 
energies have been done analytically or with the Lekner method or with the Hautmann- 
Klein method (HK)[42J. Evaluations of particle-particle interaction energies in columns 
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Ewald , refer to evaluations performed by using straightforwardly the Ewald sums for 
Yukawa potentials given by Eqs.(19), and with k x k = 16 x 16 for the number of vectors 
G in the reciprocal space and a as indicated in the Table. The term Truncation refers 
to evaluations of the energies of pairs with a direct truncation at L/2. of the Yukawa 
potentials. 

Table 4 : Average energies computed with the Monte Carlo algorithm for a Yukawa 
bilayer system. The geometrical and physical parameters of the bilayer correspond to 
run a of ref. [35J ; the coupling constant is T = 196, N = 512 is the number of particles 
and the spatial periodicity L = 28.36. f3U/N is the average total energy per particle, 
j3E int ra/N and (3Ei nter /N are respectively the average intralayer and interlayer energies 
per particle, f3au/N is an estimation of the statistical fluctuation o\j of the total energy. 
Data for the Coulomb interaction (k, = 0) are extracted from ref. [35] and data for the 
Direct Truncation have been computed with a direct truncation of the Yukawa potential 
for distance greater than L/2. 
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Systems and/or 
approximations 


K 

(inverse screening length) 


Q 

(effective charge) 


physical 
parameters 


Debye-Hiickel 

(Primitive Model 
of electrolytes) 


2 \ 9 

ek B T i 


ft 
e 

^exp(KCTy) 

e(l + KCTjj) 


Qi, Pi 
and Gij 


Colloids 
(monovalent counterions) 


k 2 = AttXbP 


Ze exp(/t<r/2) 
e(l + «<t/2) 


Z, X B , a 
and p 


Plasmas and 
Dusty Plasmas 
(a = 0) 

(a^O) 


2 47rnoe 2 ^ 1 1 ^ 
K ~ e k B [ T e + T t } 


Ze 
eo 

Zeexp(«;a) 
e (f + «a) 


Z, no, Ti, T e 
and a 



Table 1 M. Mazars, Yukawa I : Ewald sums. 
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relative accuracy 




5 ~ 10- 4 






5 ~ 10- 6 




N 


100 


1000 


10000 


100 


1000 


10000 


a* N = 10 


15.1 


1.9 


0.45 


23. 


3.1 


0.75 


ctXt — 2 


3.1 


0.7 


0.2 


5. 


1.2 


0.3 


^ = 1 


2. 


0.45 


0.1 


3.1 


0.75 


0.2 


= 0.5 


1.2 


0.3 


0.07 


2. 


0.5 


0.15 



Table 2 M. Mazars, Yukawa I : Ewald sums. 
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Configurations 




Coulomb 




Ewald 


Truncation 




Ewald 


Truncation 


and References 


methods 


Limit 


a 


K = 10~ 6 


k = 10- 6 


a 


K= 1. 


K=l. 


(2:12,2/12,212) 


















(0.5,0.5,0.) 


analytical 


-2.28472 


1.0 


-1.57705 


-1.41421 


1.0 


-0.793469 


-0.697304 


[S3 SO] 






6.0 


-2.28472 




6.0 


-1.43067 










12.0 


-2.28472 




12.0 


-1.43067 










18.0 


-2.28481 




18.0 


-1.43076 




(0.1,0.1,0.1) 


Lekner 


-5.77212 


1.0 


-5.74929 


-5.7735 


1.0 


-4.84952 


-4.85531 


001111] 


and 




6.0 


-5.77212 




6.0 


-4.87041 






HK 




12.0 


-5.77212 




12.0 


-4.87041 










18.0 


-5.77221 




18.0 


-4.8705 




(0.,0.,0.25) 


Lekner 


-3.72483 


1.0 


-3.80379 


-4.0000 


1.0 


-3.05203 


-3.1152 


@Q] 


and 




6.0 


-3.72483 




6.0 


-2.98362 






HK 




12.0 


-3.72483 




12.0 


-2.98362 










18.0 


-3.72492 




18.0 


-2.98371 




(-0.25,-0.15,-0.2) 


Lekner 


-2.82156 


1.0 


-2.73285 


-2.82843 


1.0 


-1.9636 


-1.98609 


m 


and 




6.0 


-2.82156 




6.0 


-2.04483 






HK 




12.0 


-2.82157 




12.0 


-2.04483 










18.0 


-2.82166 




18.0 


-2.04492 




(0.4,0.4,0.1) 


Lekner 


-2.28608 


1.0 


-1.81806 


-1.74078 


1.0 


-1.03453 


-0.980076 


m 






6.0 


-2.28609 




6.0 


-1.45555 










12.0 


-2.28609 




12.0 


-1.45555 










18.0 


-2.28618 




18.0 


-1.45564 





Table 3 M. Mazars, Yukawa I : Ewald sums. 
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Truncation 










0.1 


602.18 


-689.60(0) 


1291.7(8) 


0.067 


0.25 


33.910 


-191.63(7) 


225.54(7) 


0.044 


2.0 


-78.343 


-89.52(2) 


11.18(0) 


0.059 



Table 4 M. Mazars, Yukawa I : Ewald sums. 
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List of Figures 



Figure 1: Representation of D(r; a; k) and C(r; a; k) as functions of nr for several 
values of the ratio A = n/a. Function D and C are respectively defined by equations (10) 
and (B.4), the function D is used to compute the short ranged part of the energy while 
the function C is used to compute the short ranged part of the force between particles. 
As nr —>■ 0, we have D(r; a; k) ~ 2 and C(r; a; k) ~ 2, this corresponds to a non screened 
Coulomb interaction between particles. 

Figure 2: Schematic representation of the neutralizing backgrounds of three sys- 
tems : (a) monolayer with surface charge density o = —NQ/A, (b) bilayer with surface 
charge density a = —NQ/2A in each layer, and (c) slab with a volume charge density 
Po = -NQ/Ah. 

Figure 3: Representation of _E cc (Yukawa; k) for three configurations of a pair of 
Yukawa particles as functions of k. i? cc (Yukawa; k) are computed by using equation (19), 
with a = 6.0 and k x k = 16 x 16 for the number of vectors G in the reciprocal space. 
The particle 1 carries a charge Qi = +1 and is located at (0,0,0) ; the particle 2 carries a 
charge Q2 = — 1 and three positions of the particle 2 are considered in a box with periodic 
boundary conditions with L x — L y — 1.. The three positions of particle 2 are : a : (0.1, 
0.1, 0.1) ; b : (0., 0., 0.25) and c : (0.5, 0.5, 0.). Symbols refer to evaluations done with 
Ewald sums while lines represent evaluations done with a direct truncation at L/2. of the 
Yukawa potentials. For each configuration, the value of the Coulomb limit is indicated by 
a thick horizontal line and limiting values are explicitly given. These configurations have 
already been considered in some previous works [HU [34], HOI HI] (see also Table 3). 
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Figure 4: Dependence of the particle-particle interaction energy on the Ewald param- 
eter a. The configuration of the pair of particles considered are the same as those on Figure 
3, with Q\ = +1 and Qi = —1. The quantity AE cc (k, a) = E cc (n,a) — E cc (n,a = 6.) is 
computed with k x k = 16 x 16 for the number of vectors G in the reciprocal space and 
(a) k = 10~ 6 ; (b) k = 1. Inspection of the numerical values shows that, if 4. < a < 15., 
then | AE cc (n,a) |< 1(T 6 . 

Figure 5: Representation of the intralayer gn(s) and interlayer gn{s) correlation 
functions obtained for the Yukawa bilayer system with the same physical parameters as 
in Table 4. The correlation functions are defined in ref.[35]. For the Coulomb interaction, 
the correlation functions have been obtained from the ref. [35] where the computations 
have been done with Ewald, Hautmann-Klein and Lekner methods. 
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Figure 1: M. Mazars, Yukawa I : Ewald sums. 
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Figure 2: M. Mazars, Yukawa I : Ewald sums. 
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Figure 3: M. Mazars, Yukawa I : Ewald sums. 
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Figure 4: M. Mazars, Yukawa I : Ewald sums. 
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Figure 5: M. Mazars, Yukawa I : Ewald sums. 
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